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The pair-interactions Uij(r) determine the thermodynamics and linear transport properties of 
matter via the pair-distribution functions (PDFs), i.e., gij(r). Great simplicity is achieved if J7y(r) 
could be directly used to predict material properties via classical simulations, avoiding many-body 
wavefunctions. Warm dense matter (WDM) is encountered in quasi-equilibria where the electron 
temperature T e differs from the ion temperature Ti, as in laser-heated or in shock-compressed matter. 
The electron PDFs g e e(r) as perturbed by the ions are used to evaluate fully non-local exchange- 
correlation corrections to the free energy, using Hydrogen as an example. Electron-ion potentials for 
ions with a bound core are discussed with Al and Si as examples, for WDM with T e 7^ Ti, and valid for 
times shorter than the electron-ion relaxation time. In some cases the potentials develop attractive 
regions, and then become repulsive and 'Yukawa-like' for higher T e . These results clarify the origin 
of initial phonon-hardening and rapid release. Pair-potentials for shock-heated WDM show that 
phonon hardening would not occur in most such systems. Defining meaningful quasi-equilibrium 
static transport coefficients consistent with the dynamic values is addressed. There seems to be no 
meaningful 'static conductivity' obtainable by extrapolating experimental or theoretical a(uj,Ti,T e ) 
to ui — > 0, unless Ti —¥ T e as well. Illustrative calculations of quasi-static resistivities R(Ti,T e ) of 
laser-heated as well as shock-heated Aluminum and Silicon are presented using our pseudopotentials, 
pair-potentials and classical integral equations. The quasi-static resistivities display clear differences 
in their temperature evolutions, but are not the strict cj — s> limits of the dynamic values. 

PACS numbers: 52.25. Jm,52.70.La,71. 15. Mb,52.27.Gr 



I. INTRODUCTION 

Novel techniques of high-energy deposition on mat- 
ter using high-energy short-pulse lasers as well as shock 
waves enable one to produce matter in a variety of novel 
states. They demand theoretical techniques adapted to 
many physical regimes. Such warm-correlated matter 
(WCM), and warm dense mater (WDM), include hot- 
nucleating nanocrystals we well as hot-dense radiative 
(HDR) plasmas They have applications ranging 

from laser micro- machining, laser ablation [2], Coulomb 
explosions Q, inertial-confinement fusion (ICF), astro- 
physics [H and aero-space re-entry applications 5]. The 
material may be electrons and ions in a complex mixture 
of ionization states of \ow-Z and high-Z ions, e.g., H to U. 
The temperatures could be very high and yet high com- 
pressions could make the electrons degenerate or partially 
degenerate. Predictions of the thermodynamics, trans- 
port, radiative, and thermonuclear processes [l| pose ma- 
jor challenges. Furthermore, ion temperatures Tj may 
differ from the electron temperature T e , and the quasi- 
equilibrium properties, relaxation and transport have to 
be updated within the time-scales and energy-scales of 
such two-temperature equilibria Q. 

Modern electronic-structure methods provide in situ 
density-functional potentials incorporated into molecular 
dynamics (MD) simulations Q. The all-electron, all- ion 
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quantum calculations for simple systems provide 'bench- 
marks', usually at low temperatures and normal densi- 
ties. Such methods become untenable at higher-T where 
large numbers of partially occupied eigenstates have to 
be included. HDR plasmas, relevant to ICF provide a 
perspective of the problem. The complexity of WDM, 
compressions and temperatures involved call for simpli- 
fications without sacrificing accuracy. Thus molecular 
dynamics (MD), i.e., classical simulations, using various 
effective potentials have been a focus of some recent stud- 
ies [1, Q . Such methods were used in earlier work [To| , of- 
ten with phenomenological or 'chemical' models for pair- 
interactions, and in recent studies as well [l2[. Other 
authors have studied statistical potentials inclusive of 
quantum diffraction effects 13] arising from integrating 
out the electrons. A re-examination of potentials useful 
in the non-equilibrium context is needed. 

We study two types of non-equilibrium (two- 
temperature) systems, viz., (i) generated by short-pulse 
lasers (the ion subsystem remains virtually unchanged 
while T e increases during the pulse period) and (ii) me- 
chanical shocks. Here the electrons remain at T e while Tj 
and the compressions change. Our approach has been to 
construct the charge distribution n(r) around a nucleus 
immersed in the medium via the neutral-pseudoatom 
(NPA) density- functional theory (DFT) models, e.g., 
Ref. [l5[ . The n(r) is used to construct a local weak 
pseudopotential dependent on the density and tempera- 
ture of the ambient environment. The pseudopotentials 
U e i(r) provide pair-potentials Uij and pair-distribution 
functions (PDFs) gij{r). These provide the thermody- 
namics and transport properties of the system in a self- 
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consistent manner. The dependence of the f7 e j on the sys- 
tem parameters arises because the effective core-charge 
Z, the ionic-core radii r c , etc., of the ion-configuration 
depend on the medium. A type of transferable pseudopo- 
tentials is available with popular simulation software. 
However, they assume frozen-core radii r c and an ioniza- 
tion Z consistent with normal densities and T = 0. They 
are not transferable to highly compressed WDM states, 
or low electron degeneracies, as is well known and noted 
even recently by Mazevet el al. [Tlj. A more detailed 
approach that does not require defining a 'core' would 
be to use all-electron regularized pseudo-potentials based 
on norm-conserving (NC) or projector-augmented-wave 
(PAW) approaches. Such methods would be strongly 
computer-intensive and not useful for most WDM prob- 
lems. They could serve to provide a new set of bench- 
marks that are beyond Quantum-Monte Carlo calcula- 
tions, rather than serve actual WDM calculations. 

Using such NC, PAW or 'standard' pseudopotentials 
requires solving the many-center Schrodinger equation, 
as implemented in major codes [14| . The numerical sim- 
plicity needed for studying complicated WCM systems 
with many ionization states and components in differ- 
ent temperature states is lost. Hence we focus on simple 
linear- response potentials as in, e.g., Ref. [l5| . The ion 
subsystem can be treated using MD or via integral equa- 
tions like the hyper-netted-chain (HNC), or the modified- 
HNC (MHNC [16]) method to exploit spherical symme- 
try. This works even for quantum electrons at T = 
via a classical map [l7[. Hence the present study is an 
extension to non-equilibrium (two-temperature) systems 
generated by short-pulse lasers or mechanical shock. The 
physics expressed in terms of pair-potentials and PDFs 
can be directly generalized to deal with T e ^ Ti, i.e., two- 
temperature quasi-equilibria, using g su {r,T s ,T u ), where 
s, u are electron or ion subsystem labels. Such generaliza- 
tions can be formally justified in terms of the Bogolubov- 
Zubarev type of non-equilibrium theory fl8j . Within the 
time scales r e i where subsystem Hamiltonians H e and Hi 
remain invariant, we can also justify the use of DFT for 
each subsystem. 

Since laser-pulsed heating or shock-compression exper- 
iments begin at some reference density near room tem- 
perature, the pseudopotentials can be checked against 
experimental liquid-metal properties. We use such liquid- 
metal-adapted pseudopotentials, with pair potentials 
constructed from finite-T response functions incorporat- 
ing finite-T local fields consistent with the sum rules and 
finite-T exchange-correlation effects 19]. Equilibrium 
and quasi-equilibrium WDM Aluminum and Si are stud- 
ied in detail. Experimentally useful quantities accessible 
via these calculations are (a) pair-distribution functions 
and structure factors (b) thermodynamics, e.g., Hugo- 
niots, subsystem free energies, etc. (c) static and dy- 
namic liner transport coefficients, (d) energy-relaxation 
rates (e) X-ray Thompson scattering and other dynam- 
ical results. In this paper we address various aspects of 
(a)-(d) and leave (e) for a future study. 



We also ask if the electron- ion pseudopotentials J7 e i (r) 
or the pair-potentials Uu(r) could be approximated by a 
Yukawa form W SS ' , and compare them with more detailed 
potentials. 

Here we note that DFT or Car-Parincllo methods can- 
not calculate the electron-electron PDFs <? e e'M, where 
e, e' are electrons (with spin indices included as needed). 
However, these electron-electron PDFs may be calculated 
using a well-tested classical-map technique (CHNC) [13, 
|20| . inspired by finite-temperature DFT itself. This is 
used in sec. Ql] to obtain the fully nonlocal exchange- 
correlation free energy of electrons interacting directly 
with ions, rather than jellium. In sec. |TT] we use the 
electron-electron pair-distribution functions <? ee (A, r, T), 
where A is a coupling constant, to evaluate the difference 
between the exchange-correlation free energy of electrons 
in jellium and in hot hydrogen via a fully non-local eval- 
uation using a coupling constant integration of g ee (r,T) 
in hydrogen. 

Section [HT] presents the details of the pseudopotential 
and PDF calculations for ions with a core specified by the 
charge Z(Ti) and the radius r c . Results for equilibrium 
and quasi-equilibrium Al and Si WDMs are given. The 
section IIVI on transport properties addresses the mean- 
ing of a 'static' conductivity as the frequency indepen- 
dent limit of the two-temperature dynamic conductivity 
a(u>,Ti,T e ). It is argued that there may be no physical 
meaning in the 'static conductivity' obtained by extrap- 
olating a(uj,Ti,T e ) to to 0, unless T — » T e as well. 
Calculations arc presented for a possible candidate to a 
'quasi-static' resistivity of laser-heated and shock-heated 
Aluminum, showing how they differ from each other and 
from the T-dependent equilibrium resistivity. Such re- 
sults can be used to test if the physical models presented 
here provide a reasonably accurate picture of laser-pulse 
heated or shock-heated materials in the WDM regime. 



II. NON-LOCAL EXCHANGE-CORRELATION 
CALCULATIONS. 

An important application of finite-T DFT is to evalu- 
ate the exchange-correlation contribution f xc {T) per elec- 
tron to the Helmholtz free-energy / of a given electron 
distribution n(r) at the temperature T. This is usu- 
ally done in the local-density approximation (LDA) using 
the f X c(T) per unit volume of the uniform electron gas 
(UGE) as the input functional. Some authors have even 
approximated this with the T = UGE parametriza- 
tion, even though good approximations to f xc from eval- 
uations of finitc-T bubble diagrams [U [22| , the Iyetomi- 
Ichimaru (IYI) parametrization [23| as well as the f xr .(T) 
from CHNC have been available for sometime [19]. The 
IYI scheme and the CHNC f xc agree well in typical WDM 
regimes of density and T. 

In this section we present a fully non-local calculation 
of f xc at finite T for a system of electrons interacting 
with a subsystem of protons, using classical potentials 
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and pair-distribution functions as the ingredients of our 
calculation. The protons are in thermal equilibrium with 
the electrons, and have kinetic energy. However, no Born- 
Oppenheimer approximation is invoked. The ions are 
a system of classical particles and no classical map is 

needed. The proton free energy has a correlation contri- 

(v) 

bution Fc which gets automatically evaluated through 
the HNC procedure. It is a highly non-local quantity 
that cannot be evaluated in the LDA even as a first ap- 
proximation, but we do not dwell on it here. In fact, the 
total free-energy given by CHNC was used to calculate a 
Deuterium Hugoniot in jU but the nature of the purely 
electronic part f xc (r s ,T) = F xc /N, where N is the total 
number of electrons, was not examined there in detail. 
Here we high-light that aspect and concentrate on the 
electronic part of the exchange-correlation free energy 
p 

xc- 

The Hamiltonian of the system and its free energy cal- 
culation are similar to that of the electron gas, except 
that the unresponsive 'jellium' background is replaced 
by the proton subsystem. Thus 

H = H® + Hp + H ep + H ee + H pp (1) 
F = F™f +F™f + F ep + F xc + F^ ( 2 ) 

Here the superfix mf stands for 'mean-field', and all the 
terms have an obvious meaning. Since the electron and 
proton densities are identical, the electron Wigner-Seitz 
radius r s is also the proton r s . The classical map is 
needed only for the electrons. We have used the HNC 
version of the classical map [l7|, i.e., CHNC, for this 
calculation. In CHNC the quantum electrons (even at 
T=0) are treated as a Coulomb gas at a finite tempera- 
ture dependent on r s and T. The electrons interact via 
the diffraction-corrected Coulomb potential plus a Pauli 
exclusion potential. The latter accurately reproduces the 
exchange hole in the parallel-spin g e e{r) 27]. The den- 
sity (~ 0.8 g/cc) in this illustrative calculation is chosen 
so that we have a fully ionized e-p gas at a high compres- 
sion and r s — 1.5 au. The treatment of ions with bound 
states may follow the methods for Si and Al given below. 

In the panels (a) , (b) of fig [TJ we display the electron 
pair-distribution functions g ee (\,r,T) for values of the 
coupling constant A = 1, and 0.3. The Coulomb inter- 
action 1/r is replaced by A/r with A varied from to 
1. The fully non-local exchange-correlation free energy 
per electron, f xc — F xc /N, is given by an integration of 
hij( r ) — 9ij( r ) ~ 1 over A. 

fxc = l[^kj ^r[*"(r,A) + Mr,A)]. (3) 

It is noteworthy that the parallel-spin PDF, i.e., <7n(r) 
is hardly affected by the value of the coupling constant. It 
is dominated entirely by the Pauli exclusion effect. Hence 
the correlations are mostly mediated by the anti-parallel 
g(r) which changes drastically with A. The success of the 
use of an 'exact exchange' in 'optimized-effective poten- 
tials' and related methods [25[ is related to this property 
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FIG. 1: (Online color) (a)The parallel and anti-parallel elec- 
tron pair-distribution functions g ee (r) at the full value of the 
Coulomb interaction, i.e., A = 1 in a system containing inter- 
acting protons and electrons at r 3 — 1.5 au. and T = 1 eV. 
(b) The same PDFs evaluated at A =0.3, and for the same 
r 3 ,T as before, (c) The non-local exchange-correlation free 
energy f xc (solid line) in the hydrogen fluid, evaluated from 
the coupling-constant integration over the PDFs, compared 
with a similar calculation [19|] for jellium (dash-dot line) at r s 
= 1.5 and T = eV. The exchange energy f x is evaluated from 
the PDFs at A = 0. Contributions to correlations from g ep {r) 
are not included here. 



of gn(r). Spurious 'self-interaction' errors are absent in 
CHNC implementations or any methods that directly use 
a g(r). In panel (c) we present our values of f xc (H), i.e., 
the xc-free energy per electron in the presence of pro- 
tons, and the corresponding quantity, labeled f xc {J) in 
the absence of protons, appropriate for the finite-T jel- 
lium model. The exchange free energy f x is simply the 
exchange energy at A = 0, and is the same for the sys- 
tem with protons (Hydrogen) or without protons (i.e., 
jellium). Thus the correlation contribution is the many- 
body correction beyond the jellium f x at finite-T. The 
reference density matrix is diagonal in a set of plane- 
waves. The Hartree energy is zero. This definition of 
correlation energy is different to a treatment where the 
electron-proton fluid is treated as a system of hot atoms 
together with correlation corrections. There the correla- 
tion energy is usually defined with reference to a den- 
sity matrix diagonal in an atomic Hartree-Fock basis. 
The corrections to the energy beyond the Hartree-Fock 
value are deemed the correlation energy. Such a basis al- 
ready has electron-ion correlations in the wavefunctions 
which arc not plane waves. They are not included in the 
f xc used here which is purely the electronic part. In ef- 
fect, the correlation corrections in g e i(r) are differently 
partitioned in these different definitions, and hence care 
must be taken in making comparisons between atomic- 
like Hartree-Fock approaches to WDM (26j, and associ- 
ated discussions of finite-T correlation energies. 
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III. THE PSEUDOPOTENTIALS. 

Accurate pseudopotentials are available with standard 
computational packages [T3|, designed for condensed- 
matter applications where finite-T effects, multiple ion- 
izations, non-equilibrium effects etc., are uncommon. 
WCM applications require accurate but simpler poten- 
tials where classical simulations, linear response meth- 
ods, etc., could be used instead of basis-set diagonaliza- 
tions of a Hamiltonian. Here we consider (i) the simplest 
Yukawa type potentials, and (ii) more accurate pseu- 
dopotentials from DFT for a single nucleus immersed in a 
WDM medium. Thus a quantum calculation is necessary 
only for a single nucleus in a spherically symmetric envi- 
ronment modeled as a 'neutral pseudo atom' (NPA) [ll| . 
The ion density is at the ion temperature T. A DFT cal- 
culation at T e gives the electron density n(r) around the 
nucleus, the effective ionic charge Z, Kohn-Sham wave- 
functions and phase-shifts. 



1. The concept of Z. 

The core charge and the core radius r c are impor- 
tant parameters in the pseudopotentials discussed here. 
Some authors have expressed concern that Z is 'not a 
well-defined physical quantity', or that there is no quan- 
tum mechanical operator whose mean value is Z. One 
may also note that there is no quantum mechanical op- 
erator whose mean value is the temperature T. Both Z 
and T are Lagrange multipliers used to define relational 
properties. Pseudopotentials, and many other calcula- 
tional quantities have the advantage of being definable 
with respect to a convenient reference Hq, but this does 
not mean any 'arbitrariness' in the theory. The Hq se- 
lected in pseudopotential theory enables us to reduce a 
many-electron atomic problem to a few electron problem 
involving only Z electrons. 

A properly constructed Z should satisfy (i) the neu- 
trality condition n e — Zrii, where ni is the ion density. 
That is, Z is the Lagrange multiplier whose value is cho- 
sen to ensures this neutrality condition. This is readily 
generalized to a mixture of ions with different charges 
Zj (ii) the physical potential seen by a test charge in 
the plasma should tend to Z/r for large r, (iii) The Z 
should be consistent with the Friedel sum rule based on 
the phase-shifts from the ion with a charge Z and a core 
radius r c . Finally, it has to be consistent with the set 
of bound states attached to the nucleus. These issues 
are discussed in Refs. [30], [18] and provide a stringent 
set of constraints in choosing Z. We have found that 
these constraints could generally be well satisfied within 
finite-T DFT calculations that use a sufficiently large cor- 
relation sphere as the calculational volume. It has to be 
recognized that a material system, e.g., an Aluminum 
WDM, could be a mixture of many ionization states, e.g., 
Al Zj + , Zj =■■ ■ — 1, 0, 1, 2, 3 • ■ ■ etc., where all ionization 
states are integers and correspond to different descrip- 



tions of the core, with different values of r c . The concen- 
trations of each species would be such as to minimize the 
total free energy if it is an equilibrium system [l8T ] . 



2. Yukawa potentials. 

The commonly used finite-T Thomas-Fermi potentials, 
also called Yukawa potentials, are simple two-parameter 
forms. The core charge of the ions is Z and r c is set to 
zero (i.e., point ions). The field electrons are character- 
ized by a screening wavevector Ky ■ 
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We use atomic units with |e| = h = m e = 1, and the 
temperatures are in energy units. If the electron tem- 
perature is T e , the Thomas- Fermi electron- ion screening 
constant Ky is given by 



Ki 



4 

7rT R 



k 2 dkn{k){l -n(fc)} 



i(k) = 1/ [l + exp{(fc 2 /2-Ai)/T e }] 



(6) 
(7) 



Here /i is the electron chemical potential at T e for 
the electron density n e consistent with the ionization 
^(Tj,Te). In the simplest theory Z is determined via 
finite-T Thomas-Fermi theory. The structure factor con- 
sistent with this type of theory is: 

S y (k) =k 2 /(k 2 +Ky) (8) 



3. DFT-based linear-response pseudo-potentials. 

The one-electron density is the essential functional 
that determines the equilibrium properties of a quantum 
system. Hence the electron density n(r) around a nu- 
cleus of charge Z placed in a uniform electron fluid at 
the bulk density n e and temperature T e is computed. 
The computational codes for such calculations have been 
available at least since Lieberman's Inferno code [28| . 
and extended recently as in Purgatorio [23]. Most such 
codes use a single ion in a Wigner-Seitz cell of radius 
r BS = {3/(47r/9)} 1 / 3 as the computational volume when 
the ion density is p. In contrast, we use a large correla- 
tion sphere typical of the ion-correlations in the WDM. 
The single ion is surrounded by its self-consistent ion dis- 
tribution p(r) and the corresponding electron distribu- 
tion n(r), calculated self-consistently [30[. A correlation 
sphere of radius R c of 5 to 10 times r ws is used where 
9u{Rc) — ► 1. In most cases one may use the simplified 
form due to Perrot, where the ions are replaced by a uni- 
form neutralizing background except for a cavity around 
the nucleus 31| to mimic the ion-ion PDF. Further jus- 
tification of this 'neutral- pseudo-atom' (NPA) model is 
given in ref. [l5j 
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A weak local (i.e., s-wave) pseudopotential U e i(q), 
where q is the wave vector is calculated from the free- 
electron density 'pile- up' n/(r) obtained via the finite-T 
Kohn-Sham equation. This nf(q) is the density pile- up 
corrected in NPA [15j for the cavity placed around the 
nucleus to simulate the 'cavity' of the gu(r). Thus 



AU el (q) = n f (q)/xee(r s ,T e ,q) 



(9) 



Here x(r s ,T e ,g) is the electron linear- response function 
at the electron Wigner-Seitz radius r s , and temperature 
T e . A finite-T local-field correction (LFC) G q is also 
needed to go beyond the simple random-phase approxi- 
mation. That is, 



Xee(q) = X°(q)/{1-V q (l-G q )x°(q)} (10) 
G q = {l-loh}(q/k TF ) 2 (11) 
V q = 4n/q 2 , k TF = i^n/Epf/ 2 (12) 

Here x° (<?) is the non-interacting finite-T Lindhard func- 
tion. An electron effective mass m* (a 'band mass') 
is used in evaluating the Lindhard function. Thus, 
for Liquid Al near its melting point (m.pt., 0.081eV), 
to* = 0.998 gives good agreement of the calculated S(k) 
with experiment [32j . At higher temperatures to* = 1. 
The LFC G q (T e ) is taken in the LDA (i.e., the q -» 
form is used for all q), and evaluated from the ratio of 
the non-interacting and interacting finite-T compressibil- 
ities 70 and 7 respectively. Alternatively, G q with full 
q-dispersion consistent with the S(k) obtained from the 
PDF g ee (r) may be used, as discussed in Ref. [17| ■ How- 
ever, the formulation in terms of the q — > given in 
Eq. QT] is quite accurate for the systems studied here. 

Instead of numerical tables obtained from Eq. [SJ a sim- 
ple fit to the U e i(q) may be used. The fits used here are 
the Heine- Abarankov (HA) forms 33j, or a slight gener- 
alization by Dharma-wardana and Perrot (la ]. The two 
elements Al and Si chosen for this study only need a well- 
depth parameter D, and r c which is the core radius. The 
electron effective mass (i.e., band mass) to* enters into 
the response function. 

U el (r) = -ZD/r c , r < r c (13) 

= -Z/r, r>r c (14) 

U ei (q) =- ZV q M(q) (15) 

M{q) = Dsm(qr c )/qr c + (l- D)cos{qr c ) (16) 

Furthermore, it is found from test calculations that the 
scaled quantities D/(Z/r ws ), r c /r ws could be treated as 
constants for changes in r ws if Z remains unchanged. 
This enables us to explore the behaviour of these WDMs 
for compressions deviating significantly from unity, with- 
out having to do full DFT-calculations in each case. The 
use of a suitable electron-effective mass to* bypasses the 
complications of 'non-local' pseudopotentials and provide 
excellent results with just s-wave potentials. Usually m* 
can be set to unity. It also may be used as a parameter 
fine tuned to get agreement with experimental S(k) data, 



or in fitting a pseudopotential to nf(q) obtained from a 
DFT calculation, as in Eq. |H1 

The NPA approach to simple local pseudopotentials is 
reliable if the resulting U e i(q) is such that 



U ei (q)/(-ZV q ) < 1. 



A. Equilibrium pair potentials. 



(17) 



We present illustrative numerical results for Al and Si 
potentials with T e = Tj = T. The screened pair-potential 
Uu(q) is calculated from the pseudopotential via: 



Uu(q) = Z 2 (T)V q + \U ei (q)\ 2 Xee(q,r s ,T) 



(18) 



Further improvement of these formulae to incorporate 
core-polarization effects etc., via DFT are discussed in 
the appendix of Ref. [l5[ . 

Potentials for calculations with C, Si, and Ge in 
the equilibrium WDM regime have been discussed by 
Dharma-wardana and Perrot in Ref. (34[, and for Al for 
T > 5 eV in Ref. [l8[. The normal-density liquid-metal 
Al pseudopotential used here is the potential discussed in 
Ref. HI. This U ei(<l) generates the DFT n(q) under lin- 
ear response and recovers experimental S(k) data quite 
well. 

Pair potentials of other elements and at most compres- 
sions can be evaluated as in Ref. [3l|. Numerical results 
for (3Uu{r), f3 = 1/T are displayed in Fig. [5J They are 
are particularly interesting in view of the subsidiary 2kp 
peak in S(k). Solid silicon has a rather open structure 
(diamond-like) which collapses under melting into a high- 
density metallic fluid. Car-Parinello simulations, many- 
atom Stillinger- Weber type potentials (requiring many 
empirical parameters) etc., need to be highly elaborate 
to capture the subsidiary structure associated with the 
electron mass to* and the 2kp scattering processes. In 
fact, attempts were made in some of these early MD stud- 
ies to explain this 'hump' in the Si-S(k) in terms of co- 
valent bonds between Si atoms. In reality, molten Si is 
an excellent metal and the 'hump' is a metallic property 
associated with m*{k). In our view, two-center bonding 
effects do not arise unles electron densities fall to low val- 
ues (e.g., 1/64 of the normal electron density for Al), and 
the temperatures become sufficiently low. 

WDM systems generated from laser pulsing or shock 
compression usually starts from a solid-density. Thus the 
pseudopotentials and pair-potentials should successfully 
recover the strongly-correlated low-temperature regime 
near the melting point (m.pt), as well as the higher-T 
and k regimes. Our Al-Al pair potential at the m.pt 
is shown in the inset to Fig. [3J together with the DRT 
potential (based on a non-local pseudopotential) which 
is known to recover the properties of solid and liquid Al 
at low-temperatures [36| . 

It is seen that the Yukawa potentials approximately re- 
produces the envelope of the pair-potentials, while failing 
to reproduce the Friedel oscillations and energy minima. 
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FIG. 2: (Online color). (top) The Si-Si pair-potentials pUu(r), 
P = 1/T, at the melting point (T 4 = T e = T = 0.15 eV) at 
unit compression (r tus =3.07au), using an empty-core pseu- 
dopotential with core radius r c /r ws =0.3084 are displayed for 
electron effective masses m* = 1, and 1.1. The Yukawa po- 
tential's also displayed. (Bottom panel) The corresponding 
ion-ion structure factors Su(k), as well as the experimental 
data of Waseda [32| are shown. The additional structure near 
the Kohn anomaly ( 2k f) beyond the main-peak is very sen- 
sitive to the local band mass m* 



Similarly, the Yukawa structure factor S y (k) does not 
have oscillatory features. 



1. Pair-potentials for mixtures. 

Many WDM systems involve mixtures of ions. Thus 
plastic converts to a compressed WDM containing mostly 
a mixture of H, C, O ions. Even an Aluminum plasma 
could contain several stages of ionization and they have 
to be treated as a mixture. The detailed treatment of 
such systems as mixtures of average atoms, their chemical 
species-dependent potentials, embedding energies etc., as 
well as pair-potentials were discussed in Ref. [15]. The 
pair potential involving two types of ions carrying mean 
charges Z a , Z b could be written as 

U atb (q) = V q Z a Z b + \U ea (q)U eb (q)\xee(q) (19) 



Such cross-potentials would be inputs to the MHNC 
equations, or to MD simulations, for generating the cor- 
responding PDFs. 

B. Non-equilibrium pair-interactions. 

In the following we examine two-temperature quasi- 
equilibrium systems rather than full non-equilibria. Such 
systems occur during times t such that r ee < tu < t < 
T e i, after the energy is dumped in either the ion subsys- 
tem (mechanical shock) or in the electron subsystem (by 
a laser pulse). 

1. Laser-pulse generated quasi- equilibria. 

When a metal foil at the bath-temperature T b is sub- 
ject to a short-pulse laser, the electron temperature rises 
rapidly to some T e within femto-second times scales r ee , 
by equilibration via electron-electron interactions. The 
transfer of energy to the ion subsystem, arising from 
electron-ion interactions is slow, and hence the ion tem- 
perature Ti remains essentially locked to T b , while the 
electrons reach T e . Hence, processes occurring within the 
electron-ion temperature relaxation timescale r e i > t > 
Tu may be probed to provide information about quasi- 
equilibrium systems. The probes are optical pulses sam- 
pling a volume related to a space-time average over the 
pulse time and the optical depth of the material. Hence 
the 'experimental results' reported should be regarded 
as already implicitly containing some sort of interpreta- 
tional model. 

A proper description of such systems needs pair- 
potentials Uu(r,T e ,Ti) where T e ^ Ti. The procedure 
described for the equilibrium system can be simply gen- 
eralized for quasi-equilibria. In laser-heated systems, the 
pseudopotential is that of the initial state T = Ti = T e , 
and this remains unchanged under changes of T e as long 
as Z remains unchanged. Unlike in the equilibrium case, 
the bare ion-ion pair potential in a laser-heated metal is 
screened by the (hot) electrons at T e . Hence, 

U«(q,T u T B ) = {Z(T0}X + IM^)| 2 x ee (<7,r s ,T e ) 

(20) 

The ion-ion structure factors S(k) is that of the initial 
state &t T = Ti = T e for time scales shorter than r e i . 

The main part of Fig. [3] shows the evolution of the 
potentials under laser-pulsed heating where T e increases 
while the ion subsystem remains intact for time scales 
shorter than r e i. The first minimum near r/r ws ~ 1.85 
(marked ^4) weakens at first, and then deepens for T e >^ 
0.08 eV (with T, fixed at 0.03 eV). This is clearly evident 
for the case T e — 0.5 eV in the figure. This may be 
considered a microscopic signature of 'phonon hardening' 
observed in some experiments (37j . 

The highest phonon frequency ui p h supported by the 
Al-lattice at 90K is ~ 6 x 10 13 s _1 , while this de- 
creases rapidly at higher temperatures, and becomes 



7 



T.=0.03 eV 
r =2.984 au 




T 

e 


=0.03 eV 


T 

e 


=0.05 eV 


T 

e 


=0.08 eV 


T 

e 


=0.5 eV 


T 

e 


=1.0 eV 


T 

e 


=3.0 eV 


T 

e 


=10 eV 


Yu 


T =10 eV 



FIG. 3: (Online color). (Inset) The Al-Al pair-potential, 
fiiUu(r), f3i = 1/Ti, at room temperature (T e = Ti — T — 
0.03 eV) and normal density, labeled A is derived from a two- 
parameter local pseudopotential and a band mass m* — 1.02. 
It is compared with a well established non-local potential [3(| 
that accurately reproduces structural and phonon data of Alu- 
minum. The Yukawa potential, Eq. O at T = 0.03 eV is also 
displayed. The main panel displays the evolution of the po- 
tential while Ti is held fixed at 0.03 eV, and T e is increased. 
The potential becomes repulsive from T e =0.03 eV to 0.08 eV. 
For T e > 0.08 eV it deepens and becomes attractive, e.g., at 
T e = 0.5 eV, before it finally becomes repulsive. The Yukawa 
potential at T e = 10 eV is also shown 



purely imaginary if the lattice melts. The concept of 
a 'phonon' may not become meaningful if l/co p h begins 
to exceed r e j. In fact, the two-temperature pair-potential 
itself requires updating at each time step. 

When the electron temperature rises above ~ 1 eV, the 
potentials (see Fig. [3]) become repulsive and the resulting 
compression in the electron subsystem increases rapidly, 
leading to a 'Coulomb explosion'. We have assumed that 
T e is not sufficiently high to open up an ionization chan- 
nel that would add to the energy relaxation between cold 
ions and hot electrons. Such processes may be addressed 
as in Ref. HI. 



2. Shock generated quasi- equilibria 

The time evolutions of T and k in shock-compressed 
systems are quite different to those of laser-heated sam- 
ples. The Hugoniot is the locus of states that can be 
reached by a single shock wave, providing a relation con- 
necting the volume, internal energy and the compression. 
Well-developed mathematical and numerical techniques 
exist for shock studies, if the equation of state (EOS) 
is available [4fjj . The Hugoniot prevails only after many 
time steps greater than r e j. We are interested in shorter 
quasi-equilibrium time scales. Experimentally, the tem- 
perature evolution of the ion subsystem is very difficult to 
measure, and the temperature deduced from the thermal 
cmissivity of the shocked state is for the state following 
hydrodynamic interactions between the sample and the 
observation window. 

Thus material properties of the shocked state is deter- 
mined from optical measurements combined with veloc- 
ities of shock fronts and particles. The 'sample region' 
(seen by the optical probe) is the layer of material within 
approximately one optical depth behind the shock front 
(assuming that we are observing a shock wave in-flight 
inside the sample with no release wave). For a typical 
optical depth of ~ 5 nm, and shock speeds of ~ 10 4 m/s, 
the transit time T s h of the shock front through this layer 
is ~ 0.5 ps [4l|. Since r ee << T a h, there is a meaning- 
ful T e in the sample region. In any case no energy has 
been dumped into the electron subsystem for t << r e i. 
Since the phonon frequency of typical solids is of the or- 
der of several THz (i.e., Ta < T s h), one may assume that 
the ion subsystem can be characterized by a temperature 
Ti. This would lead to a single two-T quasi-equilibrium 
in the sample volume, or possibly a gradient of such 
quasi-equilibria in T e and Ti along the path behind the 
shock front. In either case, we need two-temperature 
pair-potentials along a non-equilibrium 'Hugoniot'. 

We model this system as follows. Here the usual EOS is 
a good approximation to the quasi-EOS mainly because 
the energy in the electron system is small compared to 
that of the ions. Electrons remain more or less degener- 
ate even at Ti=10 eV, when the compression reaches ~ 3. 
If k remains within a region where Z is unchanged, the 
bound electrons remain intact. However, in general the 
value of Z appropriate to the compression must be used 
and energy-relaxation by ionization need to be consid- 
ered [3!| . The compressions are calculated from Al and Si 
Hugoniots where we may neglect the small differences in 
the contribution to the internal energy from the electron 
subsystem at equilibrium or at quasi-equilibrium. The 
Al Hugoniot is based on our neutral-pseudoatom calcu- 
lations [42(| . The Si Hugoniot [43l is the L140 data based 
on the Quotidian-EOS model [44| . 

The onset of compression makes the potential take 
a Yukawa form except at very low-energy deposition. 
This pair-potential does not support phonons or 'phonon 
hardening' at the accessible short-times scales, or at 
longer time scales. In this problem the ion-electron time 
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FIG. 4: (a)The pair-potentials fiiUuir), f3 % = 1/Ti for a 
shock-heated Aluminum WDM with an initial temperature of 
T e — 0.03 eV, with the ion-subsystem temperature increasing 
to 2.5 eV. (b) The corresponding ion-ion PDFs gu(r) where 
T e remains fixed at 0.03 eV, while T% increases to 2.5 eV as 
displayed. Calculations with bridge (MHNC) and without 
bridge corrections (HNC) are given for T = 0.14 eV, showing 
that bridge corrections are negligible. The other PDFs are 
HNC solutions. The Yukawa PDF is the Boltzmann form. 



scales r e i are lengthened due to the formation of ion- 
acoustic coupled modes involving both electrons and ions 
(hence the life times of the two-T quasi-equilibria become 
longer) . 

Figure |4] shows the pair-potentials (upper panel) 
and the corresponding two-temperature pair-distribution 
functions obtained by this procedure for Al. No bridge 
corrections have been included in the g(r) shown in Fig. [4] 
except for Ti = 0.14 eV, T e = 0.03 eV, where it is shown 
that their effect is negligible. The Yukawa PDF is taken 
as the Boltzmann form exp{— WV(r)/Tj}. 

One may wonder whether potentials beyond binary- 
interactions would be necessary in a classical represen- 
tation of a many-particle system. Attempts to simu- 
late classical fluids with model potentials (e.g., Lennard- 
Jones) show that binary interactions are often inade- 
quate. However, since the potentials were made out of 



TABLE I: The temperature Ti (eV) and compression k, for 
shock-compressed aluminum and Si, for an initial normal den- 
sity and temperatures T e = T = 0.03 eV. The ionization Z 
remains essentially at 3 for Al, and 4 for Si. The electron 
temperature T e remains at the initial temperature for time 
scales less than ~ r e ». 



T(eV) 


0.14 


0.25 


0.50 


1.00 


2.50 


5.00 


K-Al 


1.4 


1.57 


1.70 


1.89 


2.26 


2.59 


K-Si 


1.57 


1.70 


1.87 


2.14 


2.65 


3.07 



Coulomb interactions in a systematic manner using DFT 
densities, we believe that such effects are not significant, 
and probably get included in the problem of modeling 
the bridge functions. However, this needs further study, 
as it could be important in low-density, low temperature 
systems, and near the critical point [45ll46|. 

IV. TWO-TEMPERATURE TRANSPORT 
COEFFICIENTS. 

The transport coefficients that are of common interest 
are associated with weak applied gradients in the elec- 
tric potential, temperature etc., so that linear transport 
coefficients like the electrical conductivity a e and heat 
conductivity ah could be defined. However, if the non- 
equilibrium system is such that the subsystem Hamilto- 
nians H e , Hi are invariant only for time scales r ee , tu, 
then the w — »■ limit of the transport coefficients, i.e., 
their static values are not defined in any strict sense. 
Thus extrapolations of experimental data to the lu — > 
have to be carefully examined to understand their phys- 
ical content. We refer to any processes with a time scale 
r such that Ta < t < r e i as quasi-static. 

Thus these experiments need to be planned, keeping in 
mind the available lifetimes for experimentation. Eval- 
uating energy-relaxation lifetimes r e i from the e-i pseu- 
dopotentials and dynamic response functions is itself an 
extremely demanding task [18j, |39( but needed in WDM 
studies. We note that if the Al-pseudopotential discussed 
here is used at unit compression, Tj = 0.1 eV, T e = 1.0 
eV, then T e j is ~ 6.02 x 10 6 electron-plasma oscillations 
long, while at Tj = 0.1, T e = 5 eV this reduces to ~ 0.72 
of the value at T e = 1 eV. These numbers are based on 
our formula for the energy-relaxation rate evaluated via 
the /-sum rule [39], without including coupled- mode ef- 
fects. 

Here we examine the theory of the dynamic two- 
temperature conductivity as a typical example. The ob- 
jective is to extract a frequency-independent (i.e., quasi- 
static) component, and present numerical calculations for 
the quasi-static part of the resistivity using the pseudopo- 
tentials and PDFs developed in the previous sections. 
This analysis puts a different perspective on the physi- 
cal content of the two-temperature Ziman formulae that 



9 




T(eV) 

FIG. 5: (Color online). Triangles: the equilibrium static resis- 
tivity of Aluminum at normal compression when T = T c = T; 
increases from the melting point. Circles: the evolution of the 
quasi-static component of the dynamic resistivity (Eq. I25p 
of WDM- Aluminum at the melting density, with the ion- 
temperature Ti held at 0.03 eV while the electron temper- 
ature T e is increased, as in laser-heated WDMs. Squares: the 
corresponding quasi-static resistivity for a shock-heated Al- 
WDM where T e remains fixed at 0.03 eV, while T increases 
to 5 eV with compression, where k increases from 1.0 to ~ 2.6 



T 




FIG. 6: (Color online). Triangles: the equilibrium static re- 
sistivity of Si at normal compression when T = T e = T in- 
creases from near the melting point. Circles: the evolution of 
the quasi-static component of the dynamic resistivity fEq, l25p 
of WDM-Si at the melting density, with the ion-temperature 
T held at 0.15 eV while the electron temperature T e is in- 
creased, as in laser-heated WDMs. Squares: the correspond- 
ing quasi-static resistivity for a shock-heated Si- WDM where 
T e remains fixed at 0.15eV, while Ti increases to 5 eV with 
compression, where k increases from 1.0 to ~ 3. 



have been used in the past for these systems. The extent 
of its physical validity has to be established by experi- 
ment. 

The response time of the system to the probe has to 
be smaller than the system relaxation time r e i , and long 
enough for the particle distributions to become station- 
ary under the probe field. If a quasi-static electric field 
E is applied to an electron distribution n(k), where k 
is a wave vector, then it is the displaced stationary dis- 
tribution n(k,E) — n(k) + S(k,E) that determines the 
electrical conductivity of the system. The probe speed 
has to be consistent with the formation of a steady- 
state displaced-electron-distribution n(k, E) in the elec- 
tron system. The probe frequency cannot be slower than 
1/Vei. Hence quantities obtained by extrapolating uj — >• 
may not have a physical meaning. 

The measured conductivity can be associated with a 



scattering time r(w) by writing 



cr(w) 



ne 2 T(ui) 



m e {l — iujT(uj)} 



(21) 



Here we have restored the electron charge and mass al- 
though we use atomic units. The conductivity o~{ui) can 
be expressed in terms of the Fourier component at the 
frequency uj of the current-current correlation function 
<< J(t),J(t r ) >>. A standard evaluation, assuming 
that the effect of the total system is a superposition of 
the effect of individual scatterers can be carried out, as 
in the equilibrium case j42}, but using the Keldysh con- 
tour for the non-equilibrium case. Here we keep in mind 
that the subsystem temperatures are T e and Ti, giving 
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the result: 



T(UJ 



q 2 dq 



2„2 



U el (q,r s ,T e )\ 2 q 



7^Im{xii(q, v, r ws ,Ti)}AN(is, w) 

xIm{l + V qX ee{q,LJ + u,r a ,T e )} (22) 
AN{v,u) = [N{{v + uj)/T e )-N{ V lT % )] (23) 



The last line contains Bose occupation factors N{{v + 
w)/T e } and N(v/Ti). Here U e i(q, r s ,T e ) is the pseudopo- 
tential presented in Eqs. The frequency-dependent 

ion response function and the electron response function 
also appear. The latter was discussed in Eq.JTU] The ion- 
response function can be constructed in a similar manner, 
or obtained from an MD simulation 48] . 

This expression for the dynamic conductivity or r(w) 
differs significantly from the Green- Kubo (GK) form used 
by a number of authors (49| since both electron-electron 
and ion-ion many-body effects are treated dynamically 
in Eq. [25] In particular, in the currently available GK 
implementations, the Green-Kubo formula is averaged 
within an ensemble of static configurations generated 
by ion molecular dynamics; the electrons are quenched 
to a Born-Oppenheimer surface to evaluate one-electron 
Kohn-Sham eigenvalues which are used in the GK for- 
mula. The latter actually calls for eigenvalues of the 
Dyson equation. The 'quenching process' destroys dy- 
namical information; the ions and electrons 'do not know' 
each other's temperatures. 

We extract a meaningful 'static' component of the 
quasi-equilibrium dynamic resistivity as this has been the 
object of recent experiments, e.g., in Ref. The phys- 
ically meaningful lowest value lul of the probe frequency 
u) is ~ l/r e i- Hence, when T e is not strongly different 
from Ti we may use uj = ljl in a Taylor expansion to 
write: 



resistivity of the form: 
p(Ti) 



R{T e ,Ti) 



e(q,r s ,T e ) 



37rn(T e ) J Q 



de— — / b ii (q,r ws ,T i ) 



de 



U ei (q,r s ,T e ) 



q dq 



e(q, r s ,T e ) 
1 + V q Xee(q,r s ,T e ) 



(25) 
(26) 



The ion- ion structure factor Su(q) is obtained from the 
ion- ion PDF via the Fourier transform of gu(r) — 1. We 
emphasize that this quasi-static component of the dy- 
namic resistivity is not the u> — > limit of the dynamic 
resistivity. 

The numerically inconvenient integration over the 
derivative of the Fermi function in Eq. [55] can be side- 
stepped by using the form of the Ziman formula give n by 
Perrot and Dharma-wardana, viz., as Eq.(31) in Ref. [15| . 

Whether the static part of the dynamic conductivity 
extracted in the above manner would be the quantity 
obtained by extrapolation of the experimental results is 
unknown. In Figs. [5] and we display the variation of 
R(Ti,T e ) for normal density Aluminum and Silicon at 
the melting point under laser heating, contrasted with 
shock heating where the variable is Tj. The compression 
remains at unity in the laser-heating and 'equilibrium' 
curves, while the compression rises in the shock- heating 
experiment. Similar results for Si are given in Fig.[B] 

It should be noted that all resistivity calculations, even 
for equilibrium systems, may differ significantly from one 
another depending on the specific (Boltzmann conductiv- 
ity or the Ziman resistivity) formula that is used. They 
themselves may differ significantly, depending whether a 
T-matrix or a pseudopotential is used [51] . If a T-matrix 
is used, corresponding modified densities of states should 
be used (here we have used the simple pseudopotential 
only). However, the trends shown in Fig. [5] and [5] clearly 
distinguish between equilibrium, laser-heated and shock- 
heated plasmas. 



DISCUSSION 



AN(v,u) = N(v/T e )-N(v/Ti)+w - K J eJ +■ • • (24) 

av 

The leading term in Eq. [24] is zero in the equilibrium 
case T e = Ti. It is seen that the first-derivative term 
gives the Ziman resistivity in the equilibrium case, where 
an uj — )• limit exists. When T e ^ Ti, the first term 
in AN is non-zero and leads to a divergence unless ui 
remains finite. Thus the concept of a static conductivity 
is non-physical unless this 'zeroth-order' correction could 
be deemed negligible. 

Replacing uj by its least allowed value uil ~ l/r e i in 
the above equation and in the pre-factor in Eq. I22[ we 
obtain a frequency-independent first-order contribution. 
This can be further reduced to a two-temperature Ziman 



We have examined the construction of effective poten- 
tials for electron-ion and ion-ion interactions that take 
account of the ambient material conditions in warm- 
correlated matter, since standard pseudopotentials are 
ill-adapted to finitc-T problems. 

The effective potentials can be used in classical MD 
simulations or modified hyper-netted-chain equations to 
obtain pair-distribution functions. The PDFs or calcu- 
lated structure factors have been used, together with the 
pseudopotentials to obtain free energies, Hugoniots, en- 
ergy relaxation and transport coefficients, even in quasi- 
equilibrium situations, with minimal computational ef- 
fort. The accuracy of the methods could be tested against 
liquid-metal WDM data. The T e ^ Ti results presented 
here have shed light on issues like phonon hardening, and 
quasi-static resistivities. 
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An alternative approach is to use statistical potentials 
generated via Slater sums 0, [53|. However, as far as 
we know, those methods have not been tested against 
experimental liquid-metal S(k) and transport date. The 
present approach based on the neutral-pseudo-atom DFT 
density is conceptually and computationally simple. The 
non-linear response of the electron system to a nucleus 
is treated by DFT, and these non-linear features are in- 
cluded in Eq. [9]where finally a linear-response pseudopo- 
tential is constructed. The weak pseudopotentials pre- 
sented here enables us to write down the pair-potential 
directly, where as the pseudopotentials provided with 
standard codes need a further solution of the Schrodinger 
two-ion problem to extract a pair potential. That ap- 
proach is indeed necessary when the condition given in 
Eq. [T7]is not satisfied. However, in the cases considered 
in this paper, the full Schrodinger two-ion problem, as 
well as the multi-ion problem (e.g., as realized in Car- 
Parinello DFT-MD calculations on C, Si, Ge molten flu- 
ids [HI], and on Al-WDM [54|) yield results in good agree- 
ment with those obtained from our single-center meth- 
ods. 

The ranges of validity of the methods discussed here 
are those of (i) the underlying DFT code used to con- 
struct the self-consistent charge density n(r) around a 
single nucleus. This is not valid if chemical bond- 
formation is possible, e.g., at low free-electron densities 
and low temperatures, (ii) the assumption that a linear- 
response pseudopotential exists, (iii) and that the rele- 
vant time scales are satisfied. High-Z materials like Au, 



W, pose special difficulties by these (or other) methods. 
The theory of equilibrium-Au WDM by these methods 
is given in Ref [521 ] . The main improvement needed is 
a relativistic-DFT calculation (55|. Once the relativistic 
charge distribution n(r) is obtained, a weak pseudopoten- 
tial may not exist. Then a two-center calculation may 
be needed, or a T-matrix has to be constructed using 
the phase shifts obtained from the DFT code and pair- 
interactions may need multiple scattering corrections. 
However, one can learn from the corresponding T = 
problem where much of the th eory is available from e.g., 
Korringa, Kohn and Rostoker [56j |. 

In laser heating or shock compression, the theory as- 
sumes that the external field has set up the and T e 
quasi-equilibrium. In particular, the laser-field is as- 
sumed to be switched off. Hence the method is indepen- 
dent of the laser intensity, polarization etc. However, if a 
strong laser field is present, then the DFT calculation for 
the charge densities and transport properties have to be 
carried out while including the self-consistent modifica- 
tion of the occupation numbers of electronic states by the 
laser field, as well as the associated dynamic screening of 
the electromagnetic field, by extending the theory given 
in Refs. [H, H3, HI] for strong electro-magnetic fields. 
Such calculations for WDM have not been attempted so 
far by anyone. 
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